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In parallel with the evolution of femtosecond and attosecond laser as well as free-electron laser 
technology, a variety of theoretical methods have been developed to describe the behavior of atoms, 
molecules, clusters, and solids under the action of those laser pulses. Here we review major ab 
initio wave-function-based numerical approaches to simulate multielectron dynamics in atoms and 
molecules driven by intense long-wavelength and/or ultrashort short-wavelength laser pulses. Di¬ 
rect solution of the time-dependent Schrodinger equation (TDSE), though its applicability is lim¬ 
ited to He, H 2 , and Li, can provide an exact description and has been greatly contributing to the 
understanding of dynamical electron-electron correlation. Multiconfiguration self-consistent-field 
(MCSCF) approach offers a flexible framework from which a variety of methods can be derived to 
treat both atoms and molecules, with possibility to systematically control the accuracy. The equa¬ 
tions of motion of configuration interaction coefficients and molecular orbitals for general MCSCF 
ansatz have recently been derived. Time-dependent extension of the A-matrix theory, originally 
developed for electron-atom collision, can realistically and accurately describe laser-driven complex 
multielectron atoms. 


I. INTRODUCTION 

Atoms and molecules, subject to visible (VlS)-to- 
midinfrared (MIR) laser pulses with an intensity typi¬ 
cally higher than 10 14 W/cm 2 , exhibit highly nonlinear 
phenomena including above-threshold ionization (ATI), 
tunneling ionization, high-harmonic generation (HHG), 
and nonsequential double ionization (NSDI) il:,[2j. HHG 
is, especially, more and more widely used as an ultra- 
short (down to attoseconds) coherent light source in the 
extreme-ultraviolet (XUV) and soft x-ray (SX) spectral 
ranges [3H53- In addition, another type of ultrashort, in¬ 
tense, coherent XUV and x-ray sources, i.e. free-electron 
lasers have emerged. These developments have triggered 
novel research activities such as ultrafast molecular prob¬ 
ing, attosecond science, and XUV nonlinear optics [6HT3] . 

Time-dependent simulations of the electronic dynam¬ 
ics in atoms and molecules still remain a challenge. 
For high-held phenomena, direct solution of the time- 
dependent Schrodinger equation (TDSE) within the 
single-active electron (SAE) approximation is widely 
used, in which only the outermost electron is explicitly 
treated, with the influence of the others expressed by 
a frozen effective potential. Naturally, however, this ap¬ 
proximation fails to treat multielectron and multichannel 
effects [DEME which are attracting increasing inter¬ 
est. Thus, various many-electron methods have recently 
been under active development. 

In this article, we give an overview of ab initio (first- 
principles) wave-function-based approaches to simulate 


multielectron dynamics in atoms and molecules driven 
by intense long-wavelength (VIS to MIR) and/or ultra- 
short short-wavelength (XUV to SX) laser pulses. Let us 
consider, within the Born-Oppenheimer (BO) (or fixed- 
nuclei) and dipole approximations, that an V-electron 
atomic or molecular system is driven by an external laser 
electric field E(f). The dynamics of the electronic system 
is governed by the time-dependent Schrodinger equation, 


i^l=H(tMt) 

where the time-dependent Hamiltonian 

H{t) = H 1 {t) + H 2 , 

is decomposed into the one-electron part, 


Hi(t) =Y^U r ifi) 

i 


and the two-electron part, 
N 


i r . * i’ 
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(1) 

( 2 ) 

( 3 ) 


( 4 ) 


for the interelectronic Coulomb interaction. h(vi,t) in 
Eq. ([3]) is given by, 

= + ( 5 ) 
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in the length gauge, with = —iVi, and, 
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in the velocity gauge, with A (t) = — J E (t)dt being the 
vector potential. 

First, in Sec. |TTJ we introduce direct solution of the 
TDSE Eq. |l]). This approach provides an exact de¬ 
scription for He, Li, and H 2 and is powerful especially 
for the investigation of multiphoton ionization by XUV 
pulses. However, its extension beyond these species is 
extremely difficult, due to exponential increase in com¬ 
putational cost. Thus, in Sec. m we discuss a general 
class of alternative methods that can handle more elec¬ 
trons, which we call multiconfiguration self-consistent- 
field (MCSCF) approaches. These are extension of those 
developed in quantum chemistry for the ground-state 
electronic state, to the dynamics involving excitation and 
ionization. Variants with different levels of flexibility are 
reviewed. Section |IV| briefly discusses alternative meth¬ 
ods that can describe multielectron atoms. These are 
time-dependent extension of the i?-matrix theory devel¬ 
oped for electron scattering from an atom and ion. 

In ab initio simulation study of multielectron dynam¬ 
ics, we usually need to (i) prepare the initial state, (ii) 
propagate the wave function in time, and (iii) read out 
physically relevant information from the wave function. 
Furthermore, since ionization is essential, it is one of the 
major issues how to treat electrons that leave the cal¬ 
culation region. Due to limitations of space, however, 
let us concentrate on the propagation of the wave func¬ 
tion in this review. We briefly note that the initial state 
can generally be obtained through either propagation in 
imaginary time or, especially in spectral methods, sepa¬ 
rate time-independent calculation of the ground state. 

Before ending the introductory section, it is worth 
mentioning other approaches that bypass explicit use 
of A-electron wave function, thus, outside the scope 
of this review. The time-dependent density functional 
theory (TDDFT) [201424] favorably scales linearly with 
N. It is, however, difficult to estimate and systemati¬ 
cally improve the accuracy of the exchange-correlation 
potential, whose form beyond the adiabatic limit is 
not yet known. Whereas an alternative called the 
time-dependent current-density functional theory [2lj 
has been proposed, only few approximations for the 
exchange-correlation vector potential are available p5] . 
Also, these methods deliver only the electron density 
or current, not the wave function, rendering the extrac¬ 
tion of physical observables difficult. Another attractive 
method that can in principle take account of correlation 
effects and extract any one- and two-particle observable 
is the time-dependent two-particle reduced density ma¬ 
trix (TD-2RDM) method [SB] [37]. A major challenge in 
this method is how to impose so-called IV-representability 
conditions, whose complete list is not known yet. Time- 
dependent quantum Monte Carlo (TDQMC) method 
[35H30] uses de Broglie-Bohm trajectories. The calcu¬ 
lation of the quantum potential or guiding waves re¬ 
quires the knowledge of the TV-electron wave function, 
thus we need some approximation. It is not obvious how 
to systematically improve such approximation, although 


Bohmian trajectories extracted from the wave function 
calculated with other methods give some insights into 
strong-field phenomena [3TH3B] . 


II. DIRECT SOLUTION OF THE TDSE 


Direct solution of the time-dependent Schrodinger 
equation (TDSE) has become possible for He, H 2 , and 
Li. A remarkable advantage of this approach is that it 
provides, in principle, an exact description. In this sec¬ 
tion, we review simulations for the se th ree species (The 
Li case is briefly mentioned in Sec. IIA). 


A. He 

TDSE simulations for He have been developed and ap¬ 
plied to study on, e.g., single- and two-photon double 
ionization [371453] as well as single ionization mmm 
ESI[56] including delay in photoemission [57]EH], and also 
on doubly excited states [59] SO] and high-field phenom¬ 
ena with a longer wavelength [BTHB5] 


1. Grid method 


Let us first describe a frequently used grid ap¬ 
proach called the time-dependent close coupling (TDCC) 
method, applied first by Taylor, Parker et al. (61]|66] and 
later by many others [37] [33] SOI SH 03] SHI EH] [53] [531 
E5US7H7D|. In the spherical coordinate system, the two- 
electron wave function T(ri,r 2 ,t) is written as, 


^(ri,r 2 ,t) 


E 

L,M,h,l 2 


p LM ( r , 




nr 2 




(7) 


where L , M denote the total orbital angular and mag¬ 
netic quantum numbers, respectively, 1 1 , Z 2 the angular 
quantum numbers of the two electrons, P^^(ri,r 2 ,t) 
the radial wave function, (i = 1,2) the combined az¬ 
imuthal and polar angles of the i-tli electron, and 

2 )= ^2 { l t rn ihTn 2 \LM)Y hmi {9. 1 )Y l2m2 {Q, 2 ) 

mi ,7712 

( 8) 

the coupled (or bipolar) spherical harmonics, with 
(?iTOiZ 2 TO 2 |LM) being the Clebsch-Gordan Coefficients. 
An orbital angular momentum eigenstate of a two- 
electron system can be specified by a combination of four 
quantum numbers either \l 1 in 1 l 2 m 2 ) (be., Yi imi Yi 2m2 ) or 
\LMI 1 I 2 ), and their mutual conversion is mediated by the 
Clebsch-Gordan Coefficients. Thus, is the explicit 
form of \LMl\l 2 ). In practice, the sums in Eq. ([7| are 
limited to a finite number of partial waves ( L , M, Zi, Z 2 )- 
If the laser pulse is linearly polarized along the 21 - 
direction, the value of M does not change throughout 
the laser-atom interaction. If we further assume that the 
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initial state has M = 0, as is the case for the ground 
state He (L = 0, S = 0), we can restrict ourselves to only 
partial waves with M = 0. Then, Eq. 0 is simplified to, 

T(n,r 2 ,f)= £ i W ri ’ r 2 ’*) (n, ,n 2 ). (9) 

r ir 2 

-L /,6 1 ,t 2 


the eigenstates as, 


$ 


LM 


(ri,r 2 )= 


n h,h,L,M 
~'ot,n i ,n 2 


x .4 


ii,Z 2 ,ni,n 2 

(ri) Fi 2 n 2 (l‘‘2 ) . 


ri 


r 2 




(13) 


Whereas Blodgett-Ford et al. m earlier used this rep¬ 
resentation to study He under a breathing mode oscillat¬ 
ing field, Parker ef a/. EB introduced it in the context 
of intense-field multiphoton ionization for the first time. 
Also, Pindzola et al. :,72i 73] applied this approach to 
electron-hydrogen scattering and later to double ioniza¬ 
tion of He and H~ [371[57] • 

By substituting Eq. ([7]) into TDSE, we obtain a set of 
coupled partial differential equations, 


*^ P hu( r i’ r 2,*) 

= ^ 5^(TM/ 1 Z 2 |H|L'M / /iZ')P^f (n,r 2 ,i), (10) 

L',M> i',Z' 

= 7)d2(n,r 2 )P^(ri,r 2 ,t) 

+ 51 V hhi[i' 2 ( r i,r2)Pi 1 (^(ri,r 2 ,t) 

,/ 0 


E 

L'M'l\ ,1: 


TTrLML'M' ( i\ t~)Lj' Al' ( >\ 

Wi^v i' (»T, r%, t)P v v (ri,r 2 ,t), 


L'M' i 


1 ’2 


( 11 ) 


called the time-dependent close coupling equations E3 
[72U73]. Here, the operators T hh , and WggfiF' 

correspond to the kinetic energy and nuclear Coulomb 
potential, the electron-electron Coulomb interaction, and 
the interaction with the time-dependent laser field, re¬ 
spectively. Their explicit forms for M = 0 can be found, 
e.g., in [551 1571 173] . Note that Tj,; 2 and v do not 
depend on M. 


It is worth mentioning that Colgan et al. have ex¬ 
tended the TDCC method described above to double and 
triple photoionization of Li mm- 


2. Spectral method 


where A denotes the antisymmetrization operator (an- 
tisymmetrizer), projecting onto either singlet or triplet 
states to guarantee the symmetry or antisymmetry of 
the spatial wave function under the exchange of identical 
particles. As basis functions Fi tU (r), one can choose 13- 
spline functions [501 . Coulomb Sturmian functions tza, 
or radial wave functions of the He + eigenstates [46] [65]. 
For the case of (singly ionized) continuum states, one 
may also use He + states for the bound electron and 13- 
splines for the continuum electron [46] . The insertion of 
Eq. (13) into the TDSE results in a coupled first-order 
ordinary differential equations for the temporal evolution 
of the coefficients C% M (t). These are, more conveniently 
after conversion to the interaction picture, integrated, 
e.g., with an explicit scheme of Runge-Kutta type [81] , 
Details are described in [331 03] |55J [77] . 


3. Example: photoionization of an excited helium atom 

As an example, let us consider single-photon single ion¬ 
ization of an excited helium atom by an attoseconcl XUV 
pulse. The excited helium has a one-electron excitation 
character to a good approximation: one of the two elec¬ 
trons is much more deeply bound than the other. Irradi¬ 
ated by an XUV pulse, say, with a photon energy of 73 
eV and a peak intensity of 10 12 W/cm 2 , predominantly 
the inner electron absorbs the photon, starts to travel 
outward, and may collide with the outer electron. We 
simulate this process with the grid-based TDCC method 
in the length gauge. 

In order to analyze the correlation-induced response 
and relaxation dynamics of the remaining ionic subsys¬ 
tem [56], we first project out all bound neutral states be¬ 
low the first ionization potential from the instantaneous 
two-electron wave function \E'(ri,r 2 ,t) and then project 
the resulting ionic part ^^(ri, r 2 , t) onto each bound 
eigenstate tp,, of He + as, 


As an alternate approach, Kamta and Starace [501 ESI 
1751177] . and some others 0TJ [330371 [751 [70] have devel¬ 
oped a spectral method, in which the time-dependent 
wave function is expanded, 

*( ri ,r 2 ,t)= 55 CPflur,), (12) 

a,L,M 

in terms of field-free eigenstates and the expansion 
coefficients C„ M (t) are propagated in time. They express 


Xi{*,t) = J ^(O^ionfor'.ijdV, (14) 

interpreted as the continuum wave function correlated 
to the ionic state ifi or the continuum wave packet in 
ionization channel i. The time-dependent population of 
channel i is given by 2 f |yj(r,t)| 2 <i 3 r. We have confirmed 
that the contribution from doubly excited (autoionizing) 
states is negligible. 

The populations of several ionic channels with the ini¬ 
tial state being \s2pAP and ls2p x 1 P are shown in Fig. [l] 
The XUV pulse is assumed to be linearly polarized along 
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FIG. 1. Time-dependent populations of several ionic states 
(channels) for He initially in the ls2p z 1 P (solid lines) and 
ls2p x 1 P (dashed lines) states. The thin dotted line represents 
the magnitude of the laser vector potential in arbitrary units. 



FIG. 2. Cylindrical coordinates (pi,zi,ipi) and (p 2 , 22 ,^ 2 ) of 
the two electrons in H 2 whose molecular axis lies on the z 
axis. 


the 2 axis and have a base-to-base pulse duration of 5 
cycles, or 280 as (see thin lines). The ls2p x 1 P state is 
composed of M = ±1, for which the radiation field op¬ 
erator is given by Eq. (8) of Ref. [57] with 

replaced by (—1) M ^ ^ ^ ' The solution is con¬ 

verged with the maximum values of L, l -\, Z 2 being 2, 5, 5, 
respectively. 


L 1 II 
0 0 0 


B. H 2 

Due to the lack of spherical symmetry, TDSE simula¬ 
tions for H 2 are more demanding than for its two-electron 
atomic counterpart He. Let us assume that the internu- 
clear axis lies on the z axis. 


1. Grid method 


In Fig. |TJ we identify two distinct time scales 156] . 
First, during the pulse (t < 280 as), instant removal of 
the Is electron leads to the population of the 2 p and 3 p 
states by shake-up. It should be noticed that the channel 
population defined above is gauge-dependent during the 
pulse , but the qualitative feature is retained even if we 
use the velocity gauge. Then, after the pulse , where the 
results are gauge-independent, population transfer from 
2 p and 3 p to other ionic states such as 2s, 3 d, and 4/ 
takes place. We refer to these delayed transitions, as 
the knock-up/knock-down processes. The motion of the 
outgoing inner electron through the cloud of the outer re¬ 
maining electron induces transitions between ionic states 
[55] . This view is further strengthened if we compare the 
results for the two different initial states. The transitions 
from 2 p and 3 p to 2s, 3d, and 4/ are clearly reduced by 
66%, 37%, and 52%, and larger population remains in 
2 p and 3 p if He is initially in the ls2p x 1 P (dashed lines) 
state compared with the case where the initial state is 
ls2p z 1 P (solid lines), since the 2 p x cloud of the outer 
electron is distributed along the x axis while the inner 
electron is ejected along the z axis, resulting in smaller 
electron-electron interaction. 


Kono et al. |82j [83] have developed a grid method with 
the molecular axis parallel to the laser polarization, by 
introducing cylindrical coordinates (p,z,ip) (Fig. [2]). If 
we define, 

, P 1 +P 2 , 1K x 

<P = <pi-v 2 , x= — 2 —’ ( 15 ) 

since the z component M (= 0,±1,±2, ••• with \ be¬ 
tween 0 and 27r) of the total angular momentum is con¬ 
served, the wave function \E '{{pj}, {zj}, {<Pj},t) (j = 1, 2) 
takes a form e lMx <&({pj}, {zj}, (f, t). Thus, the degrees 
of freedom can be reduced from six to five. Further, 
to efficiently handle the long-range nature and singu¬ 
larity of the Coulomb interaction, they have devised 
the so-called dual transformation method [521 154] , In 
this method, they introduce scaled coordinates £ and ( 
(p = f(£),z = <?(()) to replace p and z , respectively, 
and accordingly transform the wave function and the 
Hamiltonian so that the transformed wave function van¬ 
ishes at the Coulomb singular points, equidistant grids 
near the nucleus in terms of £ and ( generate small 
grid spacings in terms of p and z, and differentiation 
in the transformed Hamiltonian can be well evaluated by 
the finite difference even near the Coulomb singularity. 
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The transformed TDSE is temporally integrated with the 
alternating-direction implicit (ADI) scheme ! 8 TJ. Using 
this method, Kono et al. have simulated H 2 in an intense 
(/ ss 10 13 — 10 14 W/cm 2 ) near-infrared laser pulse to in¬ 
vestigate, e.g., formation of localized ionic states H + H _ 
|82l 185] and ionization enhanced by two-electron dynam¬ 
ics [55] . 

In order to study mainly double ionization by XUV 
pulses, Colgan, Pindzola, Robicheaux et al. [5SJ1HHHSH] 
have extended the TDCC method described in the pre¬ 
vious subsection to a molecular hydrogen, for both linear 
(along the z axis) and circular (in the xy plane) polariza¬ 
tions. In spherical coordinates the total wave function is 
expanded as, 


T(ri,r 2 ,t) = 

*£ £ 

M mi ,771.2 


1 

= WP 

P™ 1 m 2 ( r l,0l,r2,9 2 ,t) 
r i r 2 -s/sin^i v^smA^ 


i(m 1 tp 1 +m 2 ip 2 ) 

C 5 


(16) 


where Pm im2 denotes a reduced wave function, and the 
second sum is taken over toi and m 2 satisfying M = 
mi +m 2 . As is mentioned above, M is conserved for the 
case of linear polarization along the z axis. By substitut¬ 
ing Eq. (16) into the TDSE, one obtains a set of TDCC 
equations, whose explicit forms are found in [ 68 j . The 
temporal evolution of P^ lTn2 , and thus \E , (r 1 , r 2 , f), can 
be obtained through the integration of the TDCC equa¬ 
tions, e.g., by an implicit algorithm m • As variants of 
this approach, Fprre et al. [89] expand the radial and 
angular parts by B splines [81):, HD] and coupled spherical 
harmonics, respectively, as, 


^( r i,r 2 ,t) = J2 c iJk(t) B ^ r ^ Bl ^ 2 ' > yiif 2 (^ 1 ,^ 2 ), 

i,j,k ri r2 

(17) 

with k = {h,l 2 , L, M}. Schneider et al. [9T], on the 
other hand, use prolate spheroidal coordinates (£, 77 , tp), 
in which the H^" molecule is separable: 


\Jx 2 +y 2 + (z + f ) 2 + ijx 2 + y 2 + (z - § ) 2 

^ = R ’ 

(18) 

Jx 2 + y 2 + (z + §) 2 - Jx 2 + y 2 + (z - f ) 2 

V = - R -’ 

(19) 

</3 = arctan—, ( 20 ) 

x 

where R denotes the internuclear distance. They expand 
the total wave function of H 2 as, 

x ]T n mim 2 (^, m , (21) 

m 1 ,m 2 


and discretize (£,?y) by a finite-element discrete-variable 
representation. The time evolution of a reduced wave 
function H mirn2 (£ 1 , ??i, £ 2 ) V 2 , t) is calculated by use of the 
Arnoldi-Lanczos algorithm (921 . 


2. Spectral method 

Alternatively, Saenz et al. have developed a spec¬ 
tral approach [93]. The field-free eigenstates, including 
continuum, are obtained from a configuration-interaction 
calculation [94] where the Slater determinants are formed 
with one-electron Hj eigenstates expressed in terms of a 
B-spline basis in prolate spheroidal coordinates (£, 77, ip). 
They have developed their method first for H 2 with paral¬ 
lel orientation of the molecular axis to the laser polariza¬ 
tion [93] , but later extended it to an arbitrary molecular 
orientation. Using the method, Saenz et al. have studied, 
e.g., the validity and breakdown of the SAE approxima¬ 
tion for ionization and excitation yields [93[ [951 , R- and 
orientation dependence of strong-field ionization [9DD99] 
(up to ~ 10 15 W/cm intensity at 800 nm wavelength). 
Bandrauk et al. have also developed a spectral 

method with eigenstates expressed in terms of Laguerre 
and Legendre polynomials, and studied enhanced ioniza¬ 
tion of H 2 by intense ultrashort laser pulses. 


3. Vibrational degree of freedom 

Whereas the present review basically focuses on the 
electron dynamics within the fixed-nuclei approxima¬ 
tion, it is worth noting that Bachau, Martin et al. 
m\m have developed an elaborate time-dependent 
close-coupling method that treats not only the electronic 
but also the vibrational degree of freedom quantum me¬ 
chanically. In this case, one solves the seven-dimensional 
TDSE, 


^ff 0 (ri,r 2 ,l?) +V(t) - i^j W(ri,r 2 ,R,t) = 0, (22) 

where Hq denotes the field-free Hamiltonian of H 2 , and 
V ( t ) the laser-H 2 interaction Hamiltonian (these authors 
usually use the velocity gauge). Assuming negligible 
nonadiabatic couplings, i.e., Born-Oppenheimer (BO) 
approximation, we expand 4/(ri, r 2 , R, t) with fully corre¬ 
lated adiabatic BO vibronic stationary states of the form 
ip(ri, r 2 , R)x(R) with ip and x being the electronic and 
nuclear wave functions, respectively. These eigenstates 
include the bound, the resonant doubly excited, and the 
nonresonant singly ionized continuum of H 2 . The techni¬ 
cal details that make use of a B-spline basis in spherical 
coordinates are found in [9Di mm HD3| . By construc¬ 
tion, this method does not account for double ionization. 
Bachau, Martin et al. have developed their method ini¬ 
tially for linear polarization along the z direction, i.e., 
parallel to the molecular axis, but recently extended it to 
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FIG. 3. Schematic representation of the multiconfiguration 
expansion Eq. (231. Each term on the right-hand side cor¬ 
responds to a configuration $i, $ 2 , ■ • • . The Cl coefficients 
Cl,Car- are usually assumed to vary in time, whereas spa¬ 
tial molecular orbitals, represented by horizontal bars, can be 
taken as either time-independent or -dependent. 


circular polarization whose electric field is in the yz plane 
m- They have intensively studied various aspects of 
dissociative photoionization of H 2 and D 2 molecules by 
attosecond and femtosecond XUV pulses, e.g., autoion¬ 
ization mm, control by pulse duration m , electron 
localization m , and circular dichroism in molecular- 
frame photoelectron angular distributions 11041 . 


III. MULTICONFIGURATION 
SELF-CONSISTENT-FIELD (MCSCF) 
APPROACH 

Although the TDSE approach (Sec. 0 provides an ex¬ 
act theoretical framework, this method is practically un¬ 
feasible for multi electron systems beyond He, H 2 , and 
Li, especially in an intense long-wavelength laser field. In 
order to handle multi electron dynamics, time-dependent 
multiconfiguration self-consistent-field (MCSCF) meth¬ 
ods have actively been developed. The idea is to express 
the total wave function 'F(f) as a superposition of differ¬ 
ent Slater determinants or configuration state functions 

(csf) nMnni (Fi g .§ : 

P 

*(*)=£* (23) 

i 

where expansion coefficients {Ci} are called configuration 
interaction (Cl) coefficients, and bases {<!>} are the Slater 
determinants mu mg 

Xji(u) Xj2 ( r i) ••• XTv(it) 

Xli( r 2) Xj 2 ( r 2) ••• Xj«( r 2) 


Xji ( t n) Xj 2 ( r iv) ••• XTv(Uv) 

0‘i, • ■ • ,3N e I) (24) 


$l(t) = 


y/Nl 


built from N spin orbitals Xj out of 2n spin orbitals 
{4> p \p= 1,2, ••• ,n}®{a,/3} (in the spin-restricted treat¬ 
ment) with {<fip} being spatial orbital functions and a(/3) 
the up- (down-) spin eigenfunction. The summation in 
Eq. (231 with respect to configurations I runs through 


the element of a Cl space P, consisting of a given set 
of determinants. The general multiconfiguration ansatz 
Eq. (23) can represent a very wide spectrum of methods; 


whereas {Cj} are usually taken as time-dependent, they 
can also be fixed fl08j . {4> p }, and thus {$}, can be con¬ 
sidered either time-dependent or independent. The sum 
is taken over either the full-CI space Ppci composed of 
all the possible configurations I to distribute N electrons 
among the 2 n spin orbitals or any arbitrary subspace P of 
Pfci- Orthonormal spatial orbitals are often employed, 
but this choice is not mandatory. 

The time-dependent variational principle (TDVP) or 
the Dirac-Frenkel variational principle [113H115] requires 
the action integral, 


1 (^- 4 ) i*>> 


to be stationary, i.e., 


ss = - * (<a*|^> - 


(25) 


= 0, (26) 


with respect to the variations 5'!' of the total wave 
function permitted within the given multiconfiguration 
ansatz. By substituting Eq. (J23l) into Eq. ([26| and after 


laborious algebra, one can derive the equations of motion 
(EOM) for the Cl coefficients and orbital functions. 

The computational gain thanks to the use of a limited 
number of orbital functions are concisely explained in 
Ref. m ■ Let us consider a one-dimensional two-particle 
system and try to approximate its model wave function 
^(xi,X 2 ) as shown in Fig. [4] by, 


'&{x 1 ,x 2 ) ~ r. y>i(zi)<Mx 2 ), 


(27) 


where we neglect the antisymmetrization for simplicity. 
Since 'k(xi,a; 2 ) is correlated, it cannot be well described 
by the product ip{xi)<j){x 2 ) (*max = 1)- With an appro¬ 
priate number of * max , however, the expansion Eq. (27) 


can efficiently cover the major part of the wave func¬ 
tion. On the other hand, TDSE simulations as discussed 
in Sec. |Ti| explicitly treat the entire (xi,X 2 ) space, most 
of which is hardly occupied by the electrons. Thus, the 
TDSE approach is computationally much more demand¬ 
ing. Also, from Fig. [4j one can qualitatively under¬ 
stand that TD-MCSCF methods with time-dependent 
orbital functions, which move with the electron cloud, 
can efficiently represent the total wave function by a 
smaller number of orbital functions than those with time- 
independent ones. 

Before reviewing representative examples, let us em¬ 
phasize, especially to readers with a physics background, 
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FIG. 4. Correlated two-particle model wave func¬ 
tion \H(a;i, X2) expressed as multiple configurations 
'EUMxMx*) w ith single-particle orbital functions 
and (f>i(x). See text for details. 
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FIG. 5. Different representations of the identical two-electron 
wave function and orbital transformations from one to an¬ 
other. GSO: Gram-Schmidt orthogonalization, TF: Takagi’s 
factorization, UR: unitary rotation. 


that it is not very meaningful to discuss each spatial or¬ 
bital <f>p as if it were a physical entity. One can see this 
as follows by considering a two-electron system, say, He 
(Fi g .§ m\- To account for one electron ejected by a 
strong laser pulse and the other electron remaining in the 
ion, it would be reasonable to express the wave function 
as, 


Tehf = \ [ipi ( r i> i) 02 (r 2 , t) + ^ 2 (ri, t)0i(r 2 , t)] ©s, 

(28) 

with ij>i and tp 2 being nonorthogonal, where ©5 denotes a 
singlet spin function. This form is called the generalized 


valence bond (GVB) or extended Hartree-Fock (EHF) 
wave function and was used to explain the mechanism of 
NSDI 1118 ;, 1191 . The canonical orthogonalization trans¬ 
forms Eq. ( |28[ ) into an expression in terms of orthonormal 
orbitals {y>i,<p 2 }: 

2 

^N 02 =^2A i (t)\\ip i (r 1 ,t)ip i (r2,t)\\, (29) 

*=1 


where H^^jl denotes a normalized Slater determinant 
built from the direct product of a spatial orbital function 
<fii (r) and an up spin eigenfunction and the direct product 
of a spatial orbital p>j( r) and an down spin. Equation 
(29) corresponds to a two-orbital case of natural orbital 


(NO) expansion. Any arbitrary rotation, 


0 i 

02 


= u 


‘Pi 

P2 


(30) 


by a unitary matrix U yields the multiconfiguration time- 
dependent Hartree-Fock (MCTDHF) form in terms of 
orthonormal orbitals { 0 i, 0 2 }: 

2 

^mc = ^(t)||0i(r 1 ,f)0 j (r 2 ,f)||, ( 31 ) 

i,j =1 


where the matrix C with elements Cij is complex sym¬ 
metric. Inversely, Tmc can be transformed into 
through Takagi’s factorization )120] (this holds true 
for any number of orbitals). Moreover, by means of 
Gram-Schmidt orthonormalization, one can also rewrite 


Eq. (28) as, 


'Fhf+s = Hi(f)||0 , i(ri,f)0 , i(r 2 ,'f)|| 


y 2 (ri,t)0 / i(r 2 , 


(32) 


Hence, the equivalent two-electron wave function can be 
expanded in different ways Eqs. (f28l) ([29]), ([311), and 
© (Fig©. This observation emphasizes that orbital 
functions are a mathematical tool, i.e., a kind of single¬ 
particle basis functions, to construct the multielectron 
wave function, rather than a physical entity. 

An obvious extension of Eq. (29) is an expansion of the 
two-electron wave function with more than two NOs, 


*no = 'Y^A i (t)\\vi{ri,t)ipi(r 2 ,t)\\. (33) 

2=1 


The time-dependent natural orbital (TD-NO) method 
ina directly propagates A t (t) and ipi(r,t) to simulate 
laser-driven two-electron systems. We need many NOs 
to quantitatively describe correlation-induced phenom¬ 
ena such as NSDI. The time-dependent renormalized nat¬ 
ural orbital theory (TDRNOT) [121H123] propagates, in¬ 
stead of NOs, renormalized natural orbitals (RNOs) pi 
defined as, 


Pi(r,t) = \Ai(t)\(pi(r,t). (34) 

The equations of motion for RNOs can be found in |l24| . 















































A. Time-Dependent Configuration Interaction 
(TDCI) Method 

In the TDCI approach, the spatial orbital functions 
t ]) p are time-independent, and the total wave function is 
expressed by the truncated Cl expansion: 

m = * 0 c 0 {t )+e ww+E w + • • • > 

ia ijab 

(35) 

comprised of the closed-shell HF determinant <E>o built 
from the N/2 lowest orbitals, singly excited determinants 
d>“ with fa in $ 0 replaced by c/> a , 


^ E |$o>, 


(36) 


cre{t,T> 


similarly defined doubly excited determinants etc., 
with f ().) representing the up (down) spin. The expan¬ 
sion Eq. (35) is truncated at a given order. Approaches 


that include up to single excitation, double excitation, 
• • •, are called TDCI singles (TDCIS), TDCI singles and 
doubles (TDCISD), • • •, respectively. 

The orbital functions are obtained from the ground 
state Hartree-Fock method. Whereas the configurations 
$ 0 , • ■ • are fixed in time, the time-dependent Cl 

coefficients Co, Cf, • • • account for the laser-driven 
dynamics. The equations of motion for the latter can 


be derived through the substitution of Eq. (35) into 


the time-dependent variational principle. In the case of 
TDCI, which uses time-independent orbitals, the same 
EOMs can be obtained by inserting Eq. (35) directly into 


TDSE as well. The resulting EOMs are given for TDCIS 
and the length gauge by, 


d 


i-CoW = E(t).^($ 0 |f|$“)C'fW, 

a,i 

= {e a - ei )Cf(t) + E<*< l^'l*$>C$(i) 

b,j 


E(f). ($?|f|$ 0 )Co(t)- 


■E<*? 


l*|S$>C*(i) 

(38) 


is not an eigenstate of the total field-free Hamiltonian 
Hhf + H' even within the CIS ansatz. H' does not cou¬ 
ple the ground state and the singly excited states since 
(<l>o|H , | < I>“) = 0 (Brillouin’s theorem |l25j ). The detailed 
description of their implementation can be found in ma¬ 
in TDCIS, only one electron can be fully active at once; 
it is implicitly assumed that high-field phenomena are 
predominantly single-electron processes. In this sense, 
this approach can be considered as an extended ab-initio 
formulation of SAE. However, the electron can originate 
from any occupied orbital. Thus, multiple-orbital (mul¬ 
tichannel) effects are taken into account. 

In the context of high-field phenomena and attosecond 
physics, the TDCIS method was introduced |i~27l 1128] 
and implemented 126 by Santra et ai. Reference |127| 
also discusses an interesting reformulation with a concep¬ 
tual advantage in terms of a time-dependent orbital, 




Xi 


(39) 


that assembles all the single excitations from the occu¬ 
pied orbital fa. These authors later included spin-orbit 
splitting for the occupied orbitals [E2 m\- The TD¬ 
CIS has been applied to both perturbative and nonper- 
turbative multiphoton processes such as decoherence in 
attosecond photoionization m, attosecond transient- 
absorption spectroscopy - 129] , the Cooper minimum 
}131j and multielectron effects in the giant enhancement 
( 132] in HHG spectra, two-photon ionization of Ne 8+ by 
intense ultrafast x rays m, and adiabaticity and dia- 
baticity in strong-field ionization mil- 

One of the drawbacks of TDCIS is the lack of the size 
extensivity ma; the separate treatment of two subsys¬ 
tems with TDCIS, which involves up to double excita¬ 
tion in whole, is not consistent with that of the whole 
system with the same method. Also, as is usual for 


(37) approaches involving truncated expansion Eq. (35) with 


temporally fixed orbital functions, the TDCIS equations 
are not gauge invariant (see Subsection HIE). Moreover, 
its applicability is limited to systems whose initial state 
(ground state) is correctly described by the HF method. 


B. Multiconfiguration Time-Dependent 
Hartree-Fock (MCTDHF) Method 


with r = XEi r * being the dipole operator. Hhf de¬ 
notes the mean-field Hamiltonian of a one-electron nature 
that defines the HF ground state 4> 0 and orbital functions 
fa with their orbital energies £ p and that includes the 
electron-electron mean-field potential Vmf composed of 
the Coulomb and exchange operators. H' = Hi — Vmf 
accounts for the electron-electron correlation beyond the 
mean-field contribution. The second term on the right- 
hand side of Eq. (381 describes the coupling between con¬ 


figurations singly excited from different orbitals, thus can 
be viewed as interchannel interactions. It should be, how¬ 
ever, noticed that this term is present because each 4>“ 


Let us take both the Cl coefficients and orthonormal 
orbital functions as time-dependent variational degrees 
of freedom and take the sum in Eq. (23) over the full-CI 
space Pfci : 


Pfci 

= E *!(*)<?!(*)■ 


(40) 


The sum is taken over all the possible ways to distribute 
N electrons among the 2 n spin orbitals. One can use a 
more intuitive, symbolical notation, 


’I'mctdhf : ■ ■ ■ (/> n (t)} 


N 


(41) 
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FIG. 6. High-harmonic spectra from He calculated with the 
MCTDHF method for the case of 800 nm wavelength and 4 
and 8 x 10 14 W/cm 2 peak intensity. 


representing the iV-electron full-CI wave function with 
n time-dependent orbitals. In the community of high- 
field phenomena and attosecond physics, the term MCT¬ 
DHF usually refers to this most comprehensive variant, 
whereas this term has also been used for non full-CI ex¬ 
pansion in quantum chemistry. In the limit of a large 
number of n, the MCTDHF method can deliver the ex¬ 
act solutions of TDSE in principle. The EOMs for the 
Cl coefficients Ci(t) and orbital functions \fj>i (t)} ar e de¬ 
rived in Refs. [11611135 ] and discussed in Sec. | III C 
As an example, in Fig. [ 6 ] we show the high-harmonic 
spectra from He driven by a near-infrared (NIR) laser 
pulse with a wavelength of 800 nm and a peak intensity 
of 4 and 8 x 10 14 W/cm 2 , calculated with the MCTDHF 
method. The pulse is composed of single-cycle turn-on, 
constant intensity for a single cycle, and single-cycle turn¬ 
off. The simulation (in the length gauge) is sufficiently 
converged with two orbitals and the maximum angular 
momentum of 79 for each orbital. The cutoff energy E c 
calculated from the cutoff law, 

E c = I p + 3.17U p , (42) 

with I p and U p being the ionization potential and pon- 
deromotive energy, respectively, corresponds to the 65th 
and 114th order for each of the two intensities, consistent 
with the spectra in Fig. [ 6 ] It would be extremely difficult 
to obtain such spectra for NIR fundamental wavelengths 
by means of direct solution of the TDSE. 


C. Time-Dependent Complete Active-Space 
Self-Consistent Field (TD-CASSCF) Method 

The MCTDHF method, though powerful, is difficult 
to apply to large systems, since the number of the con¬ 
figurations involved (Cl dimension), and thus its compu¬ 
tational time, scales factorially with the number of elec- 


virtual Q 

{a , 6, c} 


07(d- 

?w)— 
jfe(*)- 

Mt) — £ 

M*) — l 

02 -j- 

01 — 


active A 

{t, u , v , w : x , y} 




dynamical-core 


frozen-core 


■ core C 

{hi k} 


. occupied V 

{P, 1, r , s} 


general H 


FIG. 7. Schematic explanation of the TD-CASSCF concept, 
illustrating a twelve-electron system with two frozen-core, two 
dynamical-core, and eight active orbitals. The classification 
of orbitals and the indices we use are also shown. 


trons N. Sato and Ishikawa EMI have proposed and for¬ 
mulated a more flexible method based on the CASSCF 
concept. This method introduces core (C) and active (A) 
orbital subspaces (Fig.[f|. In an intense long-wavelength 
laser field, one reasonably expects that only high-lying 
electrons are strongly driven whereas deeply bound core 
electrons remain nonionized. Accordingly, we assume 
that the nc core orbitals are doubly occupied all the 
time. On the other hand, we consider all the possible 
distributions of Na (= N — 2 nc) electrons among ua ac¬ 
tive orbitals. The CASSCF wave function can be sym¬ 
bolically expressed as, 


^CASSCF : <t>\(t)(j)\{t) ■ ■ • /> 2 c (t) 

X {^ 0 + 1 ( 1 )^ 110 + 2 ( 1 ) ' ‘ ‘ 4 > nc + riA {t)} NA - (43) 

It should be noticed that not only the active orbitals but 
also the core orbitals, though constrained to the closed- 
shell structure, vary in time, in general, responding to 
the field formed by the laser and the other electrons. 
Alternatively, it is also possible to divide the core space 
further into fixed frozen core and time-dependent dynam¬ 
ical core subspaces. The TD-CASSCF approach includes 
time-dependent Hartree-Fock (TDHF) 136 and MCT¬ 
DHF as special cases (see below). 

Hereafter, we use orbital indices {i,j,k} for core (C), 
{t,u,v,w,x,y} for active (A), {p, g,r, s} for occupied 
(V = C + A), {a, 6 , c} for virtual (Q), and {p, zz, A, 7 , 6} 
for general (% = V + Q) orbitals (Fig. [7|. For the dis¬ 
cussion of the EOMs, it is convenient to use the second 
quantization formalism, in which the Hamiltonian is ex¬ 
pressed as, 


h = J2 KK + \ £ (44) 

y.v [iis A7 

with E» = al„a va , Eft = X) CTr al a a\ T a 7T a vtT , and 

K = J ( r > P) <M r )> (45) 
































10 


% 


= dr Kir , 


<^( r l)&7 r l)#\(r2)0 7 (r 2 ) 


■ (46) 


|ri - r 2 | 

The wave function 4>i(f) for configuration I is given by, 
l$i (t)) = JJ(ol CT ) /to |uac), (47) 

ieC o- teA 

where |uac) denotes the vacuum state, and It* = {0,1} 
is an index that specifies the configuration, satisfying 
J2teA lin — NA- It should be noticed that the an¬ 
nihilation and creation operators, a and dj la , respec¬ 
tively, are time-dependent in general. 

Assuming the orthonormality of orbitals, let us also 
introduce an anti-Hermitian transformation matrix X 
whose element is defined as mm, 


XS = (<t>n\4v)- 


(48) 


Then, the temporal variation |4') of the total wave func¬ 
tion is written as the sum of the contributions from the 
evolution of the Cl coefficients and that of the orbitals: 


where, 


l*> = X>i> 6 i+ *!*>> 


x = Y / xhe?. 

[IIS 


(49) 


(50) 


The TDVP leads to the equations of motion [TOT) . llOj , 

Cl = - <$i|A-|*>, 

(51) 

{^\EjfQX - XQ^|T) = -i(A>\EjfQH - HQEjf\At), 

(52) 

where Q = 1 — P with P = |I)(I| being the projector 

onto the Cl space P. It should be noted that Eqs. (51) 


and (52) are valid for not only CASSCF but also gen¬ 


eral MCSCF wave function s with arbitrary Cl spaces 
P, including MCTDHF (Sec. Ill B) and time-dependent 


occupation-restricted multiple-active-space (Sec. HID) 
methods. 

Numerical implementation requires the transformation 
of Eq. (52) into equations of motion in terms of each or¬ 


bital function. For that purpose, we have analyzed Ejf 
and Xjf. The orbital rotation operators Ejf are catego¬ 
rized, according to the classification of the relevant or¬ 
bitals, as, 

{k} = {e),EI El, El El, e;, El) . (53) 


By noting the projector Q in Eq. (521, we can classify 


these into the following two, for the case of TD-CASSCF, 
MCTDHF, and TDHF: 


1. Redundant (E), E* u , E(f). Both E£\V) and E£\V) 
lie inside P or vanish. In this case, Eq. (52) reduces 


to a trivial identity 0 = 0, and the corresponding 
Xjf may be matrix elements of an arbitrary one- 
electron anti-Hermitian operator 9(f) j llBj . 


Xjf = (4>im |<M, e\t) = -0(t). 


(54) 


This redundancy originates from the fact that the 
unitary rotation inside the core (C), active (A), or 
virtual ( Q ) orbitals does not alter the space cov¬ 
ered by the given MC ansatz Eq. (23). This ob¬ 


servation reemphasizes that molecular orbitals are 
mathematical instruments rather than physical en¬ 
tities. The simplest choice is 9(t ) = 0, thus Xjf = 0. 

2. Non-redundant uncoupled (Ej, E\, Ef,, El). At 
least one of Ejf\^) and Ejjf^) do not vanish, and 
Ejf |T) and Ejf |SP) lie, if non-vanishing, outside P. 
This type of rotations do not contribute to the 
Cl equations, Eq. (51), through the second term. 


Nevertheless, they contribute to the orbital EOMs 
Eq. (52), which reduces to a simpler expression 


nns|. 

(T| 


El El 


\*)Xj = -i(*l 


Ejf, H 


|4/). (55) 


Explicit formulas for Xj, Xj, Xf and Xjj are given 
in Eqs. (33)-(36) of Ref. [IU9j, with Rjf = iXjf. 

The orbital EOMs explicitly in terms of orbital func¬ 
tions, required for numerical implementation, are given 
by, 

= —iQF\(j)i) + E l<M A 7 > ( core ) ( 56 ) 

V 

I f>t) = -iQFtlcft) + E l<M A ?> (active) (57) 


for core and active orbitals, respectively, where, 

Q = i-E^^I’ 

V 

is the orbital projector onto the Q space, and, 

m) = f i^+YDuG^i), 

tu 

Ft\<h) = M> + E w'Mp:: (. d - 1 )*: 

uvwx 

f\<t> p ) = h\<t> p ) + 2 Y / G j j \<t> P ), 


= WM - l/2W p \<f> u ), 


with, 


WJ((r) = J dr 1 


,<^( r ')<M r') 


(58) 

(59) 

(60) 
( 61 ) 
(62) 

(63) 
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being the general mean-field operator, and D l u = 
<'F A |£ t “|* A > and P*” = (^aI^TI^a) being one- and 
two-electron reduced density matrix (RDM) elements, re¬ 
spectively, defined within the active space. The second 
terms on the right-hand sides of Eqs. (|56|) and (57) repre¬ 


sents, aside from the redundancy within the same orbital 
subspace, the coupling between the core and active or¬ 
bitals. We need not explicitly work with virtual orbitals 

thanks to Q m- 

The wave function Eq. (23) with configurations Eq. 


(471 can be rewritten as, 


|*)=$c|^a), (64) 

with the active part, 

p 

|4/a> = £ |I)Ci, (65) 

I 


where, 


with Ec = 2 fj ■ Whereas both options, just shifting 
the origin of energy, are equivalent, the former is numer¬ 
ically more stable [109 : . 

The TD-C ASSCF method is gauge invariant (see Sub¬ 
section HIE) and size extensive. Accuracy can be sys¬ 


tematically controlled between TDHF and MCTDHF. 
The account of correlation can also be flexibly controlled. 
Moreover, correlation (beyond the HF mean-field poten¬ 
tial) in the ground state can be taken into account, unlike 
in the case of the TDCI method. The number of the con¬ 
figurations involved (Cl dimension) scales factorially with 
the number of active electrons Na, significantly reduced 
from that in the MCTDHF approach, which scales facto¬ 
rially with the total number of electrons N. The detailed 
discussion on the computational cost of the TD-CASSCF 
method can be found in Ref. [l09j . 


D. Time-Dependent Occupation-Restricted 
Multiple-Active-Space (TD-ORMAS) Method 


I 1 ) = II n ( a L) 7t,T \vac), (66) 

a teA 


and the core generator, 

( 67 ) 

iec 


As is already mentioned above, MCTDHF and TDHF are 
special cases of the TD-CASSCF method. The former 
corresponds to nc = 0 , for which <f>c becomes the iden¬ 
tity operator, and the latter to ua = 0 and |Ta) = | vac). 
Therefore, the equations of motion for the Cl coefficients 
and orbitals presented in this subsection are valid also 
for MCTDHF and TDHF. 


The separation of the core wave function in Eq. (64) al¬ 
lows us to reduce the equations of motion for the Cl coef¬ 
ficients, Eq. (511 to those in terms of the active-part wave 
function 4^ / |109i llTOj . i.e., effectively an iVA-electron 
problem 


Ci = ~i(I\H a ~ £ A i|* a) - <I|*|*a>, (68) 

h a = flK + * £ ’ ( 69 ) 

tu tuvw 


where 1 is the identity operator, E A = (4'aIHaI’I'a), and 

/* = <^I-/K>, (70) 

g%u = {<k\-wZ\M- (7i) 


Here, without loss of generality, we have adopted a partic¬ 
ular phase choice (H/1 \R) = 0. For another, more common 
choice i(4'|T) = Eq. ( 68 ) is replaced by, 


Ci = -i{I\HA + E c i\*A) - (im’J'A), (72) 


The classification into core and active orbital subspaces 
introduced in the TD-CASSCF reduces the CI dimension 
compared with that in the MCTDHF. Nevertheless, the 
computational cost increases exponentially with Na, thus 
large-active-space calculations will become prohibitive. 
Also, we cannot take a large core subspace if we are to 
properly describe inner-shell excitation and ionization by 
short wave-length pulses. We can further decrease the CI 
dimension by resorting to non-complete CI expansion. 

The ORMAS model [ llOi 457] offers a highly flexible 
framework in this direction. In this model, we further 
subdivide the active orbital space A into a given number, 
G, of subgroups, 


G 


■A = 

s=l 


(73) 


where, 

A g = {$\$\ ■ ■ ;</>%} (1 <9< G), (74) 

with n A = n g- The minimum and maximum num¬ 

bers of electrons that can be assigned to each subgroup 
are specified as, 

Nf n < N g < iV“ ax (1 < g < G), (75) 

with N a = Eg=i N g. 

This ansatz involves a wide variety of CI expansions, 
some of whose examples are shown in Fig. [8] 

1) CASSCF: G = 1, K nn = N^ x = N A [Fig. § (a)]. 

2) HF+X: G = 2, N A - L < N ± < N A , 0 < N 2 < L 
for a given L, m = N A /2 [Fig. [j| (b)]. The CI 
space is composed of the HF reference space (Mi) 
and up to L-fold excitation from A\ to A 2 - Let us 
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(a) CASSCF (b) HF+X (c) RAS 


FIG. 8. Schematic examples of the ORMAS-CI space for six 
active electrons and twelve active orbitals. The vertical ar¬ 
rows represent electrons in the HF configurations, to be dis¬ 
tributed following the respective ORMAS restriction. The 
curved upward arrows represent the excitation from one to 
the other subgroup. 


denote the corresponding method, explicitly con¬ 
sidered by Miyagi and Madsen [13814140] , as HF+X 
for brevity, with X = S, SD, SDT, etc. It should be 
noted that this expansion is distinct from TDCI 
methods, in that the orbital functions are time- 
dependent. 

3) RAS: G = 3, and the maximum number of holes 
Mho\e in and the maximum number of electrons 
M e i ec in A 3 are specified, whereas IV2 is uncon¬ 
strained, i.e., 


Na — Mhole < Ni < N\, 0<N 3 < Melee- 


(76) 


This model is called restricted active space (RAS) 
model proposed by Olsen et al mu. 


Note that Miyagi and Madsen FT!) use a term “TD- 
RASSCF” for their method that uses the second type 
of Cl spaces above but that this method is actually not 
based on the RAS scheme of Ref. H2U, introduced much 
earlier. For consistency with the terminology widely 
used in the stationary electronic structure theory and 
to promote collaboration between atomic physicists and 
quantum chemists, in this review, we refer to the latter 
method as RAS, which includes the former as a special 


case. 


The equations (51), (52), (56), (57), and (68) of mo¬ 


tion themselves hold also for the TD-ORMAS approach. 
For the case of TD-CASSCF, including MCTDHF and 
TDHF, the orbital rotation operators are classified 
as either redundant or non-redundant coupled, as dis¬ 
cussed in the previous subsection. General TD-ORMAS 
methods, on the other hand, require a third category as 
follows: 

1. Redundant. {Ej,E£} as well as active intra-group 
rotations {E*; </> t , ()>„ G A s } belong to this. Equa¬ 
tion (52) is reduced to Eq. (54). 


2. Non-redundant uncoupled. {E\, E\, E%} belong 

to this, as in the case of TD-CASSCF. Equation 


(52) is simplified to Eq. (55). 


3. Non-redundant coupled. Either |$) or £7|\I/) 
lies across P and Q. Active inter-group rotations 
{££;</>* G Ag^u G Ag<,g ± g'} belong to this, 
in general. These contribute to both the Cl and 
orbital EOMs. We need to work directly with 
Eq. ( [52] ) . The coupled equations to be solved for 
G A g ,(f)u G A g ',g ^ g'} are given in 
Eq. (57) of Ref. [Tffil . 

Numerical implementation of the TD-ORMAS method 
is described in Ref. 1101 . It is worth noting that the or¬ 
bital EOMs Eqs. (pTand ([57]), and Cl EOMs Eq. p|, 


derived by Sato and Ishikawa in Ref. m are valid even 
for general MCSCF wave functions with arbitrary Cl 
spaces P if we treat all the active-active rotations E (j as 
non-redundant coupled and solve Eq. (57) of Ref. illOj 
for X*. The TD-ORMAS method is, however, advanta¬ 
geous from the viewpoint of the efficiency and stability of 
time propagation, since it can limit non-redundant cou¬ 
pled rotations only to active inter-group rotations. Equa¬ 
tions of motion for general MCSCF wave functions have 
recently been derived independently also by Haxton and 
McCurdy IT32| . 


E. Remark on the Gauge Dependence 


The TDSE Eq. ([I]) can be represented in either the 
length Eq. ([5]) or the velocity gauge Eq. §. The gauge 
principle states that all physical observables are gauge 
invariant [143] . The wave functions \Hl(£) and 'Fy(t) ex¬ 
pressed in the length and velocity gauges, respectively, 
are related to each other through the gauge transforma¬ 
tion, 

*v(t)=W(i)tf L (t), (77) 


where the transformation operator, 


U(t) = exp 


N 

—iA(t) ■ r i 

i —1 


(78) 


is unitary. One can easily verify this, by substituting 
Eq. ( [77] ) into the TDSE with Eq. (|6| and finding that 'Ll 
indeed satisfies the TDSE with Eq. ([5]). 

Let us denote the orbital functions calculated with 
a given multiconfiguration ansatz Eq. ( |23[ ) within the 
length gauge by {^(r)}. Equation (77) is fulfilled if one 
constructs the wave function Ty (f) of the same ansatz 
with the Cl coefficients unchanged using the orbital func¬ 
tions { 4 >p (r)} defined by, 


4>p (r) = exp [~iA(t) ■ r] <%(r). (79) 

Since this requires the orbital functions to be time- 
dependent, TD-MCSCF methods with time-independent 
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orbital functions such as TDCI discussed in Subsec, fill Al 
are, in general, not gauge invariant, i.e., the results ob¬ 
tained within the length gauge are not equal to those 
within the velocity gauge. This is because 'X'v(i) does not 
necessarily belong to the subspace of the Hilbert space 
spanned by { < f>i}, in which ’I'L(i) is optimized. 

It should be noticed that the length- and velocity- 
gauge Hamiltonians Hp (t) with Eq. (|5[) and Hy(t) with 
Eq. §, respectively, are related by [ 1431 , 

H y = + i . (80) 
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Then, since the gauge-transformation operator U{t) is 
unitary [Eq . (|78| l, one can see that the formulas for the 
TDVP Eq. (|26|) in the two representations are equivalent. 
This guarantees that the wave function transformed via 
Eq. (791 from the wave function satisfying the length- 
gauge TDVP satisfies the velocity-gauge TDVP. There¬ 
fore, TD-MCSCF methods with orbital functions varying 
in time are gauge invariant in general [109[ 11101 139 


F. Time-Dependent Coupled-Cluster method 

A somewhat different direction of particular interest is 
time-dependent extensions |1441I145| of the so-called cou¬ 
pled cluster method, though they have not been applied 
to laser-induced ionization dynamics yet, and further the¬ 
oretical sophistications appear to be required. Whereas 
coupled cluster is often abbreviated as CC, we do not use 
it in this review to avoid confusion with close coupling in 
Sec. nn 

Let us first assume that orbital functions are fixed 
in time, as implemented by Huber and Klamroth 1144] 
to study laser-driven excitation of small molecules. We 
write the total wave function as, 


FIG. 9. Partitioning of the radial coordinate of the ejected 
electron in the Zi-matrix approach. 


In order to simulate dynamics involving ionization 
without the need of a huge basis set, it would be advan¬ 
tageous to allow orbital functions to vary in time (cf. dif¬ 
ference between TDCI and HF+X). Kvaal has developed 
orbital adaptive time-dependent coupled-cluster method 
1145] , though Ref. jll5j assumes a time-independent 
Hamiltonian, thus no external field. This method is 
based on bivariational (rather than variational) principle; 
|4') and (4''| are independently variated, and constructed 
based on different sets of molecular orbitals {</j p } and 
{<p p }, respectively, that are not individually orthogonal 
but satisfy biorthogonality {(p p \ip q ) = S pq . The EOMs 
for the variational degrees of freedom, i.e., {<p P }, {^ P }, 
and their corresponding cluster amplitudes are derived 
from the principle that the action integral is stationary 
under all possible variations, as detailed in Ref. [1451 . It 
seems that the initial ground state cannot be obtained 
by imaginary time relaxation. 


IV. A-MATRIX APPROACH 


|4/(t))=e f |4/ 0 ), (81) 

where |4'o) denotes a reference Slater determinant, taken 
to be the HF ground state, and the cluster operator T is 
defined as, 

T = T i + ■■■■> ( 82 ) 

ia ijab 


which is truncated at a given term in practice. The time- 
dependent cluster amplitudes rf, , ■ ■ ■ correspond to 
C“, C'°f , ■ ■ ■ in Eq. (351 in the first order, whereas 
Eq. (81) involves excitation of any arbitrary order even if 
we truncate Eq. (82), e.g., at the second term. As a con¬ 
sequence, the exponential form recovers size extensivity 
which the TDCI method (Sec. III A| ) lacks. If we include 
up to the second term in Eq. (82) (coupled cluster sin¬ 
gles and doubles) and substitute the cluster expansion 
Eq. (81) into the TDSE, we obtain, after some algebra, 
the EOMs for the amplitudes r“(t) and r“ fc (t), whose 
explicit forms are given in Ref. nn. 


In this section we briefly mention approaches based 
on the A-matrix theory originally developed for electron- 
atom and electron-ion collision processes. The length 
gauge is preferred in these approaches [ 146] . 

In the time-dependent R-matrix (TDRM) approach, 
originally proposed by Burke and Burke [1471 and further 
developed by van der Hart, Lysaght, et al. [T381 1149] . 
single ionization at most of an (N + l)-electron atom 
is considered, and the time-dependent wave function 
^(ri, • • • , r/v+i, t) is expanded as [147] , 

T( r i,--- ,rjv+i,£) 

= M^$i(ri, • • • , r N ; f2jv+iCTjv+i)?'jv+ 1 V , i( 7 'JV+i, *), 

(83) 

where denote time-independent channel functions 
formed by coupling the residual A-electron state with 
the angular and spin eigenfunctions of the ejected elec¬ 
tron, and ipi the time-dependent reduced radial functions 
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that describe the motion of the ejected electron in the i- 
th channel. 

By introducing an equally spaced temporal mesh, 


t q = qAt, q = 0,1,2, 


(84) 


with At being the time step, one can express the tem¬ 
poral evolution of the wave function 4/(t) from t = t q to 

lq+1 &S, 


[H(t q+1/2 )-EMt q+1 ) = Q(t q ), 
with t q+ 1/2 = t q + At/2, E = 2*Af _1 and, 


(85) 


©(^) — j+1/2) + E]^(t q ), 


( 86 ) 


where terms of 0(At 3 ) are neglected. 


In order to solve the implicit step Eq. (85), we parti¬ 


tion the radial coordinate r/v+i of the ejected electron 
into an internal and external regions with the boundary 
radius do (Fig. |9|. In the internal region, exchange and 
electron-electron correlation effects between the ejected 
and the remaining electrons are important, whereas in 
the external region the ejected electron moves in the local 
long-range multipole potential of the residual iV-electron 
system and the laser field. The external region is further 
partitioned into p subregions with the boundary radii 
(do <)di < • • • < d p (Fig. pi. The outermost boundary 
dp is set so large that i/>i(a p /t) = 0. The matrix equations 
that couple ipi(a s ,t ) and (/i'(a s ,l) (s = 0, • • • ,p), where 
0' = between neighboring boundaries are obtained 

and solved as described in Refs. )147l 1149) . Whereas 
the TDRM theory was originally formulated in terms of 
the i?-matrix basis functions [H9J, also a mixed finite- 
difference and i?-matrix basis set technique has recently 
been developed, in which finite-difference methods are 
adopted for the external region 11501 . 

Guan et al. (15l| l!52j have developed a time- 
dependent extension of the R-spline R-matrix (BSR) 
method |153H155j , for multiphoton single and double ion¬ 
ization of atoms. The time-dependent wave function T 
is expanded as, 

T(ri, • • • ,r N+1 ,t) = ^2C q (t)$ q { ri, - ,rjv+i), (87) 


in terms of a known set of time-independent (N + 1)- 
electron states and expansion coefficients C q (t) de¬ 
scribing the dynamics. are formed from appropriately 
symmetrized products of atomic orbitals that are not nec¬ 
essarily orthogonal, by means of the BSR method. The 
temporal evolution of C q {t) is obtained by the Arnoldi- 
Lanczos approach 


V. CONCLUSION 


Today, thanks to elaboration of quantum chemistry 
methods, not only specialized theoreticians but also ex¬ 
perimentalists can routinely calculate ground-state prop¬ 
erties of large systems containing tens to hundreds of 
electrons. In marked contrast, approaches for the dy¬ 
namics in a strong, time-dependent external field is still 
in the early stage of development. Major difficulty stems 
from spatially extended continuum electrons and dynam¬ 
ical electronic correlation more prominent than the static 
one in the ground state. In this paper, we have reviewed, 
if not exhaustively, major promising activities to tackle 
this grand challenge. Whereas its extension beyond three 
electrons would not be realistic, full TDSE simulation 
(Sec. 0 has been a powerful method to investigate the 
electronic correlation in photo excitation and ionization, 
and will continue to be a stringent benchmark for any 
new approach. Furthe r de velopment and sophi stica tion 
in both MCSCF (Sec. Ill) and R-matrix (Sec. IVI ap¬ 


proaches may someday realize routine calculations for 
large complex atoms and molecules. In the present re¬ 
view, we have concentrated ourselves on the electron dy¬ 
namics within the fixed-nuclei approximation except for 


in Sec. IIB 3 Comparison with experiments will involve 
nuclear motion. The existence of different dissociation 
and Coulomb explosion channels may imply that nuclear 
dynamics also should be, at least partially, treated quan¬ 
tum mechanically, for which another endeavor such as in 
Ref. 1159] will be required. Vast unexplored terrain lies 
before us. 
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